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Abstract 

An efficient algorithm is developed to construct disconncctivity graphs by a random walk 
over basins of attraction. This algorithm can detect a large number of local minima, find energy 
barriers between them, and estimate local thermal averages over each basin of attraction. It 
is applied to the SK spin glass Hamiltonian where existing methods have difficulties even for a 
moderate number of spins. Finite-size results are used to make predictions in the thermodynamic 
limit that match theoretical approximations and recent findings on the free energy landscapes 
of SK spin glasses. 

PACS numbers: 05.10.Ln, 75.10.Nr, 02.70.Rr 

Disconnectivity graphs (DGs) [HE], widely used for representing energy landscapes, summarize 
local minima and energy barriers of an energy function into a tree. The DG of a continuous energy 
surface can be constructed by computational approaches that search for local minima and saddles 
based on the gradient and the Hessian matrix [3]. However, such approaches cannot be applied to 
Ising Hamiltonians defined on discrete spins. For example, the Hamiltonian of the SK spin glass [1] 
with zero external magnetic field, which is the focus of this paper, is 

H(s) = - ^JjjSjSj, (1) 

i<j 

where s = (s±, . . . , sn), Si € {±1}, is a vector of N spins and Jjj is the interaction between Sj and 
Sj. Many studies have been conducted, such as in [2H2]> to characterize the free energy landscape 
of the SK spin glass by investigating solution structures of the TAP free energy equations [ID] . 
See [IT] for a review and more references. These studies reply heavily on specific assumptions for 
the distribution of the disorder J = {Jij}- On the other hand, computational approaches have 
been developed to construct DGs for spin glass Hamiltonians jl2H18] . From DGs one may extract 
microscopic information to characterize free energy landscapes. In principle, these approaches can 
be applied given any possible distribution of the disorder, but, unfortunately, they are feasible only 
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for small-scale systems with less than or around 30 spins, due to the computationally expensive 
nature of DG construction. 

The purpose of this paper is to develop an efficient algorithm that is able to construct DGs 
containing hundreds of local minima for spin systems with N on the order of 100 and possibly larger. 
The algorithm is motivated by the broad success of the Wang-Landau (WL) algorithm |19|I20| which 
produces a random walk in energy space. To build a DG, we aim to generate a random walk over 
the basins of attraction of local minima. Suppose that the Hamiltonian H (s) has K local minima, 
v\, . . . , vk- The basin of attraction of v^, denoted by D^, is the set of configurations which will be 
sent to Vk by steepest descent that recursively flips the single spin giving the maximum decrease 
in H(s). If a random walk can be produced over all basins of attraction, D%, . . . , Dk, not only do 
we have all the local minima but also may estimate local thermal averages over every Dk- Such 
estimation on basins of attraction is a key to the utility of the inherent structure approach |21[I22[ 
and the superposition approach [23] • Furthermore, frequent transitions between basins must occur 
during the walk. We say two configurations x and y are neighbors, denoted by x y, if they differ 
by only one spin. As each local move is a single-spin flip, cross-basin moves can be used to find the 
barrier between two basins defined as B^e = mm p£ -p ke max sgp H(s), where Vki is the collection of 
all paths between and vg (k ^ £) in the configuration space. For a spin system, 

B M = mm{H(x) V H{y) : x G D k ,y G D e ,x <-> y}, (2) 

where H(x) V H(y) = max[H(x),H(y)]. Thus, keeping track of cross-basin moves, we may obtain 
a rough estimate of B^e which can be refined by a ridge descent algorithm to be introduced later. 
With local minima and barriers detected constructing the DG of H (s) is trivial. 

Since the number of minima increases exponentially with N for SK spin glasses, we follow 
the practical convention to construct DGs with K lowest local minima for a big K. For the sake 
of understanding, we first describe the algorithm in the context that K local minima of H(s) 
have already been detected. These local minima are used to partition the space into K basins, 
D\,... ,Dk, and their complement Dq. As energy of SK spin glasses is continuous, a ladder of 
energies, no < u\ < ■ ■ ■ < ul = oo, where uq is a lower bound of H(s), are employed to partition 
the energy space into L intervals. Then, our goal is to generate a random walk over all (nonempty) 
subregions, D^j = {s € : H(s) G [uj-i, Uj)}, k = 0, . . . , K and j = 1, . . . , L, where two indices, 
the basin index k and the energy index j, are used for space partition. The desired random walk 
can be implemented by a generalized WL (GWL) algorithm [241125] . where energy on a subregion 
-D/y is not a constant. Let O^j denote the (unnormalized) statistical weight of D^j in the Boltzmann 
distribution, i.e., 6kj oc Y^D k e~^ HlyS \ where j3 is the inverse temperature. A flat histogram over all 
Dfrj can be produced if the probability of visiting s G D^j is proportional to e~^ H<yS ^ /O^j- Since O^j 
is unknown we set 0^- = 1 at the first iteration of the walk. At iteration t, let 6^ be the estimate 
of 8kj and xt and y be the configurations before and after a randomly chosen spin is flipped. A 
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steepest descent operation is applied on y to find its basin index, in which the energy change of a 
single-spin flip can be computed efficiently by utilizing the additive structure in ([I]). If y is not in 
the basin of any of the K minima, then y G Dq. In general, if x± G D k j and y G the Metropolis 
ratio from x± to y is 

r(x t -+y) = rnin{l,e^)-« 0®/6®} . (3) 

Each time a subregion Dfcj is visited, the weight #Jy will be updated to 0^*- = / with a 
modification factor / > 1. Following the WL algorithm, / is reduced to v7 when the flatness of 
the histogram becomes acceptable (maximal fluctuation < 25%) since the last reduction of /. If 
the energy ladder is dense enough such that the energy in [uj-x,uj) is approximately a constant, 
then the local density of states Q k j oc 0jye* u * _1 , where Qkj is the number of configurations in the 
basin D k with energy Uj—\. In this scenario, local thermal averages over a basin of any temperature 
can be obtained via estimated 9k j, similar to the calculations in [18H20] . 

Suppose the random walk has been simulated for n iterations. Let (xt,Xf) be a pair of 
configurations simulated at two consecutive iterations, i.e., \t — t'\ = 1. For any two basins D k 
and Di we keep track of the configuration pair 

(a,6)fe£ = arg min {H(x t ) V H(xf) : x t G D k ,x t > G D e }, 

(x t ,x t ,) 

for 1 < t, t' < n. At the last iteration, (a, b)i-e is the pair that minimizes (|2|) among all cross-basin 
neighbors generated by moves between the two basins, which provides a rough estimate of B k e- A 
ridge descent algorithm is developed to refine the estimate. Let (ao, bo) = (a, b)ke such that ao G 
and 6o £ D^. For i = 1,2,..., find iteratively 

a t = arg min{ H(a) : a G Ngb(6 t _i) n D k }, 

a 

b t = argmin{i?(6) : b G Ngb(a t ) fl D e }, 

b 

until b%—\ = bt, where Ngb(s) is the set of all the neighbors of s. This iterative algorithm moves 
(ao,&o) downhill along the ridge separating the two basins. For every pair of k and £, the barrier 
B^i will be estimated by H(at) V H(bt) at the final iteration of the ridge descent. 

Next we discuss how to identify K local minima in the burn-in period of the random walk. 
A collection of local minima, V, is dynamically accumulated using the same GWL update ([3]) in 
the burn-in period, but with a constant modification factor / = e. Initially, V is empty. In every 
iteration, a local minimum is located by steepest descent starting from the proposed configuration 
y. Then, V is updated to include the lowest K minima identified so far or all of them if there are 
less than K minima identified. With In / = 1, ln^*- simply records the number of visits to D k j. 
This makes it easy to update these weights when the local minima in V are updated. As pointed 
out in many previous studies, the WL update with a big / enables the walk to reach all subregions 
very quickly, which is the key to detecting sufficient low-energy minima. 
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The proposed algorithm may be tested on the SK spin glass model ([I]), where the Jij are 
independent Gaussian random variables with mean and variance 1/N. Hereafter, only rescaled 
energy (energy per spin) will be used. We constructed DGs for small-scale systems with N = 25 
for which exact results can be obtained via enumeration as well as larger N where enumeration is 
impossible. For each N the algorithm was applied to 100 independent samples of J = {Jij} with 
P = 1/T C = 1. A rough energy range of interest of this model is [—0.8,-0.3]. Accordingly, the 
energy space was partitioned into L = 10 intervals with Uj = —0.8 + j Au for j = 0, . . . , 9, where 
Au = 0.1 for N < 60 and Am = 0.05 for N > 70. 

For N = 25, our algorithm was applied to each sample with a total of 1 x 10 7 MC sweeps. We 
chose K = 500 which turned out to be greater than the total number of minima for all the samples, 
ranging from 56 to 310. Compared to results from enumeration, the constructed DGs were highly 
accurate. Our algorithm did not miss a single minimum for any sample. Recall that due to the 
use of steepest descent our algorithm will not produce any false minima. The average absolute 
energy difference between estimated and exact barriers was 1.4 x 10~ 7 , which was extremely small 
relative to the energy range of the model. This demonstrates that the algorithm indeed accurately 
recovered most energy barriers. In fact, our algorithm recovered exactly all the barriers for 99 out of 
the 100 samples. Finding barriers is a difficult job especially for discrete Hamiltonians. The result 
here highlights the advantage of simulating a random walk over basins of attraction in building 
DGs. The average acceptance rate for the MC moves was ~30%. Thus, by exploring only 9% of all 
configurations our algorithm was able to construct DGs almost identical to those by enumeration. 

We applied our algorithm to N = 40, 50, . . . , 100, aiming at constructing DGs for the lowest 
K = 500 minima. Each run consisted of 5 x 10 8 MC sweeps. The acceptance rate for MC sweeps 
was > 15% for each N, averaging over the samples. At the final iteration, In / decreased below 
10~ 6 for most of the samples with N < 60 and was on the order of 10 -5 to 10 -3 for N > 70 
(Table [I]). These results suggest that our algorithm well explored identified basins and made 
frequent transitions between them, which is sufficient for constructing DGs although estimation of 
the weights Okj may not be very accurate for large N. Figure [1] shows the DG constructed for a 
sample with N = 100. One sees two almost identical subtrees, each containing a ground state and a 
few groups of local minima, joining at the highest detected barrier. The identical structure between 
the two subtrees, due to the fact that H{—s) = H(s) ([I]), gives a validation of the DG. However, 
the algorithm did not recover the energy landscape for those missing high-energy minima. This 
limitation is inevitable due to the exponential increase in the complexity of SK spin glasses. To 
quantify the statistical error of a constructed DG, we applied independently our algorithm to this 
sample ten times. Remarkably, all the identified local minima and at least 95% of the estimated 
barriers were exactly identical between any two runs. Furthermore, we systematically compared the 
two subtrees of a DG to measure the accuracy of our algorithm. For all N, the two subtrees of every 
sample contained identical sets of local minima, up to reversal of all the spins, and substantially 
overlapping sets of barriers (Table [T]), which demonstrates the reliability of the constructed DGs. 
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Table 1: Convergence and accuracy for N > 40. Note: — log(ln/) is the negative logarithm (base 
10) of the median, over 100 samples, of In / at the end of simulation; 7] is the percentage of barriers 
that are identical between the two subtrees of a DG, averaging over 100 samples. 
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Figure 1: A constructed DG for the SK spin glass with iV = 100. A terminal node (leaf) on the 
tree represents a local minimum and an internal node (branch point) represents an energy barrier, 
with energy levels given by the vertical axis. The statistical error between independent runs is very 
small (see text for more discussion). 

Define the barrier height h of a minimum as the energy difference between the minimum and 
its nearest barrier (its parent on the tree) . We grouped minima according to their energy u relative 
to the global minimum u* and studied the relation between (h) and iV for each group, where {X) 
denotes the average of X over samples of the disorder. We analyzed five groups of minima with 
(u — u*) £ [0.01(2; — l),0.01z) for z = 1,... ,5. In each of these energy intervals, our algorithm 
detected more than 1000 minima over the 100 samples for N = 100. A power law, (h) = cN x , 
was fitted with extremely high R 2 (> 0.98) for each group, where R is the correlation coefficient 
between ln(/i) and In A. The high consistency across different A serves as a confirmation for the 
accuracy of this result. Figure E^a) shows the fitted power laws for three groups (z = 1, 3, 5), from 
which we see the three lines are almost parallel to each other and that (h) clearly decreases with 
the increase of the energy of a minimum. The estimated A for z = 1, ... ,5 were —1.70 ± 0.06, 
-1.63 ± 0.07, -1.54 ± 0.09, -1.47 ± 0.08 and -1.53 ± 0.09, respectively. These powers were not 
significantly different especially between neighboring groups. This result implies that the barrier 
height of a minimum vanishes as N — > 00 and the rate of decay is comparable among minima with 
different energy. Thus, all minima become marginally stable as N — > 00, which is consistent with 
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Figure 2: Results for N > 40. (a) Average barrier height (h) and its fitted power laws for three 
groups of minima with energy ranges given by the intervals in the legend. Log scale is used for both 
axes, (b) (P(|g|)) for T = 0.4. Solid lines with symbols are results for N = 40,70 and 100 from 
this study. Dashed and dotted lines are, respectively, the result for iV = 128 obtained by Monte 
Carlo simulation and the result of Parisi's prediction as N — > oo, both from [26J. 

the recent finding that each minimum and its nearby saddle on the free energy surface of the SK 
spin glass coalesce in the thermodynamic limit [HE]. 

A constructed DG may also provide qualitative understanding about pure states. When the 
temperature T is low, it is reasonable to approximate a pure state a by a local minimum and 
determine its statistical weight by w a tx e - Nu a/ T j where u a is the energy of the local minimum. 
Then, one may find the probability distribution, P{q) = Yl a 1 w a w 1 5{q ai — q), for the overlap q ai 
between two pure states, a and 7, and its average, (P(q)), over samples of the disorder. From the 
constructed DGs we approximated (P(q)) for T = 0.4. Since the distribution is symmetric around 
q = 0, we plot (-P(|<ji)) in Figure E(b). Our approximation is compared against direct Monte Carlo 
simulation for N = 128 and Parisi's prediction as N — > 00 under the same temperature [26]. With 
the increase of N, the overlap distribution from our approximation becomes closer to the expected 
shape and the location of the peak for N = 100 is in good agreement with Parisi's prediction. This 
shows the utility of DGs in characterizing the key features (e.g., the order parameter) of spin glasses 
for low temperature. However, when the temperature is high, say close to T c , the statistical weight 
of the missing high-energy minima will be larger and this approximation is likely to underestimate 
-P(|<zD for small \q\. 

Computational approaches that combine local optimization and Monte Carlo sampling, such 
as this work, have been developed for global optimization with applications to protein and peptide 
models [27H29] . These existing methods were not designed to construct DGs for Ising spin models 
and are different in nature from this work. In addition, the basin-sampling approach [30J employs 
the WL algorithm to construct the total energy density of states, which shares some common 
features with the present work. Although we have focused on the SK spin glass with Gaussian 
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interactions, it should be noted that our algorithm is applicable to other possible choices of the 
disorder J and many other spin systems. For a continuous system, our method can be employed to 
find local minima with a suitable local optimization algorithm and provide rough energy barriers 
which may be refined with alternative geometry optimization methods [3]. 

I thank the two referees for their helpful suggestions and comments. This work was supported 
by NSF grant DMS-0805491. 
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